home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 12 / Mac Magazin and MacEasy Magazine CD - Issue 12.iso / Sharewarebibliothek / Anwendungen / Wissenschaft & Technik / Yorick / yorick11-ppc folder / include / rkutta.i < prev    next >
Text File  |  1995-04-03  |  8KB  |  279 lines
  1. /*
  2.    RKUTTA.I
  3.    4th order Runge-Kutta integrator (rkutta)
  4.    Also Bulirsch-Stoer integrator (bstoer)
  5.    Both after routines in Numerical Recipes by Press et.al.
  6.  
  7.    $Id$
  8.  */
  9. /*    Copyright (c) 1994.  The Regents of the University of California.
  10.                     All rights reserved.  */
  11.  
  12. /* ------------------------------------------------------------------------ */
  13.  
  14. func rkutta(derivative, y0,x0, x1,epsilon, dx0)
  15. /* DOCUMENT y1= rkutta(derivative, y0,x0, x1,epsilon, dx0)
  16.      integrates dydx= DERIVATIVE(y,x) beginning at (X0,Y0) and
  17.      going to X1 with fractional error EPSILON.  The result is
  18.      the value of y at X1.  DX0 will be used as the initial guess
  19.      for a step size.
  20.  
  21.      If the external variable rk_nstore is non-zero, rk_y and rk_x
  22.      will contain up to rk_nstore intermediate values.
  23.  
  24.      The external variable rk_maxits (default 10000) is the
  25.      maximum number of steps rkutta will take.  The variable
  26.      rk_minstep (default 0.0) is the minimum step size.  The
  27.      variable rk_maxstep (default 1.e35) is the maximum step
  28.      size, which you may need if you are storing intermediate
  29.      values (particularly with bstoer).
  30.  
  31.      If a function rk_yscale(y,dydx,x,dx) exists, it is used
  32.      to compute an appropriate yscale to give the EPSILON error
  33.      criterion meaning.  Otherwise, yscale is taken to be:
  34.         abs(y)+abs(dydx*dx)+1.e-30
  35.  
  36.      Based on odeint from Numerical Recipes (Press, et.al.).
  37.  
  38.    SEE ALSO: bstoer, rk_nstore, rk_maxits, rk_minstep, rk_maxstep,
  39.              rk_ngood, rk_nbad, rkdumb, rk4
  40.  */
  41. {
  42.   extern rk_x, rk_y, rk_ngood, rk_nbad;
  43.   rk_ngood= rk_nbad= 0;
  44.   if (rk_nstore) {
  45.     if (rk_nstore<2) rk_nstore= 2;
  46.     rk_x= array(double(x0), rk_nstore);
  47.     rk_y= array(double(y0), rk_nstore);
  48.     s= 1;
  49.   }
  50.  
  51.   dxsign= sign(x1-x0);
  52.   dx= double(abs(dx0)*dxsign);
  53.   x= double(x0);
  54.   y= double(y0);
  55.   for (n=1 ; n<=rk_maxits ; ++n) {
  56.     dydx= derivative(y, x);
  57.     if (!is_void(rk_yscale)) yscale= rk_yscale(y,dydx,x,dx);
  58.     else yscale= abs(y)+abs(dydx*dx)+1.e-30;
  59.     if (rk_nstore) s= rk_store(y,x,s);
  60.     if (abs(dx) > rk_maxstep) dx= dxsign*rk_maxstep;
  61.     if (dxsign*(x+dx-x1) > 0.0) dx= x1-x;
  62.     local dxdid, dxnxt;
  63.     y= rkqc(y,dydx, x,dx, derivative, epsilon,yscale, dxdid,dxnxt);
  64.     x+= dxdid;
  65.     if (dxdid == dx) ++rk_ngood;
  66.     else ++rk_nbad;
  67.     if (dxsign*(x-x1) >= 0.0) {
  68.       if (rk_nstore) {
  69.     s= rk_store(y,x,s);
  70.     rk_y= rk_y(..,1:s);
  71.     rk_x= rk_x(1:s);
  72.       }
  73.       return y;
  74.     }
  75.     if (abs(dxnxt) < abs(rk_minstep))
  76.       error, "required step less than rk_minstep";
  77.     dx= dxnxt;
  78.   }
  79.  
  80.   if (rk_nstore) {
  81.     s= rk_store(y,x,s);
  82.     rk_y= rk_y(..,1:s);
  83.     rk_x= rk_x(1:s);
  84.   }
  85.   error, "exceeded rk_maxits iterations";
  86. }
  87.  
  88. local rk_nstore, rk_maxits, rk_minstep, rk_maxstep, rk_ngood, rk_nbad;
  89. /* DOCUMENT rk_nstore, rk_maxits, rk_good, rk_bad
  90.      rk_nstore      maximum number of y values rkutta is to store
  91.         after calling rkutta, rk_y and rk_x will contain stored values
  92.      rk_maxits      maximum number of steps for rkutta (default 10000)
  93.      rk_minstep     minimum step size for rkutta (default 0.0)
  94.      rk_maxstep     maximum step size for rkutta (default 1.e35)
  95.      rk_ngood       number of good steps taken by rkutta
  96.      rk_nbad        number of failed (but repaired) steps taken by rkutta
  97.  */
  98. rk_maxits= 10000;
  99. rk_minstep= 0.0;
  100. rk_maxstep= 1.0e35;
  101. rk_nstore= 0;
  102.  
  103. func rk_store(y,x,s)
  104. {
  105.   ++s;
  106.   if (s > rk_nstore) {
  107.     y2= rk_y(..,1:0:2);
  108.     x2= rk_x(1:0:2);
  109.     s= numberof(x2);
  110.     rk_y(..,1:s)= y2;
  111.     rk_x(1:s)= x2;
  112.     ++s;
  113.   }
  114.   rk_y(..,s)= y;
  115.   rk_x(s)= x;
  116.   return s;
  117. }
  118.  
  119. func rkdumb(derivative, y0,x0, x1,nsteps, nosave=)
  120. /* DOCUMENT y_of_x= rkdumb(derivative, y0,x0, x1,nsteps)
  121.      integrates dydx= DERIVATIVE(y,x) beginning at (X0,Y0) and
  122.      going to X1 in NSTEPS 4th order Runge-Kutta steps.  The
  123.      result is dimsof(Y0)-by-(NSTEPS+1) values of y at the points
  124.      span(X0, X1, NSTEPS+1).
  125.      If the nosave= keyword is non-zero, the returned value will
  126.      simply be the final y value.
  127.  */
  128. {
  129.   dx= (x1-x0)/nsteps;
  130.   ++nsteps;
  131.   if (!nosave) y= array(double(y0), nsteps);
  132.   for (i=2 ; i<=nsteps ; ++i) {
  133.     y0= rk4(y0,derivative(y0,x0), x0,dx, derivative);
  134.     x0+= dx;
  135.     if (!nosave) y(..,i)= y0;
  136.   }
  137.   return nosave? y0 : y;
  138. }
  139.  
  140. func rkqc(y,dydx, x,dx, derivative, epsilon,yscale, &dxdid,&dxnxt)
  141. {
  142.   x0= x;
  143.   y0= y;
  144.  
  145.   for (;;) {
  146.     x= x0+dx;
  147.     if (x==x0) error, "integration step crash";
  148.     /* first take two half steps... */
  149.     dx2= 0.5*dx;
  150.     x2= x0+dx2;
  151.     y2= rk4(y0,dydx, x0,dx2, derivative);
  152.     y2= rk4(y2,derivative(y2,x2), x2,dx2, derivative);
  153.     /* ...then compare with one full step... */
  154.     y1= rk4(y0,dydx, x0,dx, derivative);
  155.     /* ...to estimate error */
  156.     y1= y2-y1;
  157.     err= max(abs(y1/yscale))/epsilon;
  158.     if (err <= 1.0) {
  159.       dxdid= dx;
  160.       dxnxt= (err>6.e-4)? 0.9*dx*err^-0.20 : 4.*dx;
  161.       break;
  162.     }
  163.     dx*= 0.9*err^-0.25;
  164.   }
  165.   return y2 + y1/15.;
  166. }
  167.  
  168. func rk4(y,dydx, x,dx, derivative)
  169. /* DOCUMENT y_at_x_plus_dx= rk4(y,dydx, x,dx, derivative)
  170.      takes a single 4th order Runge-Kutta step from X to X+DX.
  171.      DERIVATIVE(y,x) is a function returning dydx; the input DYDX
  172.      is DERIVATIVE(y,x) at the input (X,Y).  This fourth evaluation
  173.      of DERIVATIVE must be performed by the caller of rk4.
  174.  */
  175. {
  176.   dx2= 0.5*dx;
  177.   x2= x+dx2;
  178.   dydxp= derivative(y+dydx*dx2, x2);   /* slope at 1st trial midpoint */
  179.   dydxm= derivative(y+dydxp*dx2, x2);  /* slope at 2nd trial midpoint */
  180.   dydxp+= dydxm;
  181.   dydxm= derivative(y+dydxm*dx, x+dx); /* slope at trial endpoint */
  182.   return y + (dydx+dydxp+dydxp+dydxm)*(dx/6.0);
  183. }
  184.  
  185. /* ------------------------------------------------------------------------ */
  186.  
  187. func bstoer(derivative, y0,x0, x1,epsilon, dx0)
  188. /* DOCUMENT y1= bstoer(derivative, y0,x0, x1,epsilon, dx0)
  189.      Bulirsch-Stoer integrator, otherwise identical to rkutta routine.
  190.      All of the options for rkutta (rk_nstore, etc.) work here as well.
  191.  
  192.    SEE ALSO: rkutta, rk_nstore, rk_maxits, rk_minstep, rk_maxstep,
  193.              rk_ngood, rk_nbad
  194.  */
  195. {
  196.   extern _rzextr_x, _rzextr_d;
  197.   rkqc= bsstep;
  198.   _rzextr_x= array(0.0, numberof(_bs_nseq));
  199.   _rzextr_d= array(double(y0), 7);
  200.   return rkutta(derivative, y0,x0, x1,epsilon, dx0);
  201. }
  202.  
  203. func bsstep(y,dydx, x,dx, derivative, epsilon,yscale, &dxdid,&dxnxt)
  204. {
  205.   x0= x;
  206.   y0= y;
  207.  
  208.   for (;;) {
  209.     for (i=1 ; i<=numberof(_bs_nseq) ; ++i) {
  210.       n= _bs_nseq(i);
  211.       y= rzextr(i, (dx/n)^2, mod_midpt(y0,dydx, x0,dx, derivative, n),
  212.         yerr, 7);
  213.       err= max(abs(yerr/yscale))/epsilon;
  214.       if (err < 1.0) {
  215.     dxdid= dx;
  216.     if (i==7) dxnxt= 0.95*dx;
  217.     else if (i==6) dxnxt= 1.2*dx;
  218.     else dxnxt= (16./*_bs_nseq(6)*/*dx)/n;
  219.     return y;
  220.       }
  221.     }
  222.     /* step failed, claimed to be unusual */
  223.     dx*= 0.0625;  /* related to numberof(_bs_nseq) */
  224.     if (x+dx == x) error, "integration step crash";
  225.   }
  226. }
  227.  
  228. _bs_nseq= [2, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96];
  229.  
  230. func mod_midpt(y,dydx, x,dx, derivative, nstep)
  231. {
  232.   dx/= nstep;
  233.   ym= y;
  234.   y+= dydx*dx;
  235.   x+= dx;
  236.   dydx= derivative(y,x);
  237.   dx2= 2.*dx;
  238.   for (--nstep ; nstep ; --nstep) {
  239.     swap= ym + dydx*dx2;
  240.     ym= y;
  241.     y= swap;
  242.     x+= dx;
  243.     dydx= derivative(y,x);
  244.   }
  245.   return 0.5*(ym+y+dydx*dx);
  246. }
  247.  
  248. func rzextr(iest, xest, yest, &yerr, nuse)
  249. {
  250.   _rzextr_x(iest)= xest;
  251.   if (iest==1) {
  252.     _rzextr_d(..,1)= yest;
  253.     yerr= yest;
  254.     return yest;
  255.   } else {
  256.     m1= ((iest<nuse)? iest : nuse) - 1;
  257.     fx= _rzextr_x(iest-1 : iest-m1 : -1)/xest;  /* no more than nuse-1 */
  258.     yy= yest;
  259.     v= _rzextr_d(..,1);
  260.     _rzextr_d(..,1)= c= yy;
  261.     for (k=1 ; k<=m1 ; ++k) {
  262.       b1= fx(k)*v;
  263.       b= b1-c;
  264.       ok= double(b!=0.0);
  265.       bad= 1.0-ok;
  266.       b= ((c-v)/(b+bad))*ok;
  267.       ddy= c*b + bad*v;
  268.       c= b1*b + bad*c;
  269.       if (k!=m1) v= _rzextr_d(..,k+1);
  270.       _rzextr_d(..,k+1)= ddy;
  271.       yy+= ddy;
  272.     }
  273.     yerr= ddy;
  274.     return yy;
  275.   }
  276. }
  277.  
  278. /* ------------------------------------------------------------------------ */
  279.